
if (!file.exists("output")) {
  dir.create("output")
}

# load data
source("HNC_metadata_tidy.R")
source("data_handling.R")
source("risk_factor_regressions.R")

sbsCOSMIC <- read.csv("Input_signature_attributions/Input_SBS_COSMIC_attributions.csv", stringsAsFactors = T, row.names = 1) %>% tibble::rownames_to_column("donor_id")
dbsCOSMIC <- read.csv("Input_signature_attributions/Input_DBS_COSMIC_attributions.csv", stringsAsFactors = T, row.names = 1) %>% tibble::rownames_to_column("donor_id")
idCOSMIC <- read.csv("Input_signature_attributions/Input_ID_COSMIC_attributions.csv", stringsAsFactors = T, row.names = 1) %>% tibble::rownames_to_column("donor_id")

additional_df <- list(sbsCOSMIC, dbsCOSMIC, idCOSMIC)

cosmic <- additional_df[[1]]
for (n in 2:length(additional_df)) {
  cosmic <- cosmic %>% left_join(additional_df[[n]], by = "donor_id")
}

# merge with cosmic sigantures and tidy epi variables
data <- data %>%
  left_join(cosmic, by = "donor_id") %>%
  mutate(
    subsite = factor(subsite, levels = c("Larynx", "OC", "OPC", "Hypopharynx")),
    age_group = factor(age_group, levels = c("55-65", "0-45", "45-55", "65-75", "75+")),
    age = as.numeric(scale(age_diag)),
    alcqty = scale(ifelse(alcohol_ever == "Non-drinker", 0, alcqty))
    )

# Select signatures present in >n% of samples
# Dichotomize the signatures (Presence (signatures>0)="1",absence (signatures==0)="0")
# If signature is present in more than 75% cases, dichotomize by the median
data <- dichotomizeSigs(metadata = data, sigsn = cosmic, sigcutoff = 2)

# Dichotomized signatures to analyze
signatures = c("SBS16_cat","DBS4_cat","ID11_cat")
parameters = "tob_alc"
confounders = c("age","sex","region","subsite","alcqty")

data <- data %>% filter(!is.na(alcqty))

regression <- LRmodel(data, signatures, parameters, confounders, pval = 1, make_plot = F)

write.csv(regression[["regression_param"]], "output/Supplementary_Table_13.csv", row.names = F)

